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1 Introduction 

Recently two of the present authors Q reported on a method for computing safe bounds for the 
value of the Titchmarsh-Weyl m function associated with the differential expression 

My=-(-(py')' + qy) (1.1) 

w 

defined over [a, oo), — oo < a, where p,q,w are real-valued functions which satisfy p" 1 ,^, w £ 
Lj oc [a, oo) and w(x) > a.e. In the case that w = 1 Weyl Q showed that the differential equation 

My = Ay, A G C+ U C_ (1. 2) 

has at least one solution that belongs to the set 

POO 

L^(a, oo) = {/ : / w \ f \ 2 dx < oo}. 

J a 

The proof of this result introduced the Titchmarsh-Weyl m function to the mathematics literature. 

It was however Titchmarsh who investigated the properties of m(A) as an analytic function of A 

and established the connection between the location of its poles, of necessity on the real line, and 

the eigenvalues of the differential equation ( [1. 2j ) . 

The analytic form of the m function is determined by the form of the I? solutions of ( |1. 2j ) and 

it is perhaps not suprising that there are few examples of m known in closed form. For example 

if a = and p = w = 1 then, when q = ±£ 2 , the m function is known as a rational function of 

gamma functions, while if q = ±x it may be written in terms of Bessel functions. For a detailed 

discussion of the m function, together with examples, see ||. 

*34E05, 34E20, 34L05, ODE's, Spectral theory, Titchmarsh-Weyl m- function, verified computation 



In p] the central role of the m function in the spectral theory of Ql. 1[ ) is established. It is 
shown that the behaviour of the m function near the real line classifies all points of the real line as 
belonging to the point spectrum, continuous spectrum, point-continuous spectrum or the resolvent 
set of the self-adjoint operator generated in L^(a, oo) by Ql. 1|) together with initial conditions. A 
further use for the m function is in determining the best constant in Everitt's HELP inequality. 
See || for further details. 

In view of the central importance of the m function and also in view of the difficulties in obtaining 
its values analytically, much effort has been expended in devising computational algorithms to 
estimate its value. These are reported on in a series of papers which include ||,@,||. However 
these papers only contain estimates for the value of m and do not seek to address the question of 
absolute bounds on the error in the computations. 

In |l| we reported on an algorithm to compute rigorous bounds for the m function. This algo- 
rithm worked well for examples q = x a , a > 2 or a = 1, but was shown to be computationally 
inefficient for examples q = —x a , a = 1,2, and the interval based algorithm needed in the compu- 
tation in |l| did not cover the case q = ±x a , < a < 2, a 7^ 1. The purpose of this paper is to 
present two new algorithms which overcome these difficulties and enable the m function now to be 
enclosed for a much wider class of problems. 

In section 2 we review the relevant extracts from the theory of the Titchmarsh-Weyl m function 
that are needed to develop our algorithm. Section 3 is devoted to recalling certain asymptotic results 
that are central to our method as well as presenting an overview of an interval based algorithm 
that is fundamental to the implementation of our m computation. Section 4 contains the results 
of the numerical experiments that we have performed, while section 5 deals with the extension of 
the algorithm to overcome the problems encountered with q = ±x a , < a < 2, a ^ 1. 

2 Titchmarsh-Weyl limit-point, limit-circle classification 

In the classical limit-point, limit-circle theory of ( |l. 2|) it is shown that, starting from a pair of 
solutions 6, 4> of Ql. 2| ) which for strictly complex A satisfy 

6(a,\) = { P 6')(a,\) = l 

0(a,A) = -l (W>')(a,A) = 0, (2. 1) 

there exists a complex-valued function m(A) such that 

M; A) = 0(; A) + m(X)4>(; A) G L 2 w (a, 00). (2. 2) 



When, up to constant multiples, there is precisely one solution of Ql. 2| ) in L^(a, oo), we say ( p.. 1| ) 
or ( |1. 2| ) is limit-point at infinity. If all the solutions of 2j ) are in -C/^(a, oo), we say that ( |1. 1| ) or 
2j ) is limit-circle at infinity. Further, the limit point, limit-circle classification is determined by 
p, q, w and is independent of the strictly complex parameter A. In this paper we shall be exclusively 
concerned with the limit-point case. The m function is a Nevanlinna function, mapping the upper 
(lower) half-plane to itself, and as such has any singularities confined to the real line. From (|2. 1| ) 
and ( |2. 2| ) it follows that 

where ip is any (non-zero) constant multiple of V'o- The result ( |2. 3| ) is the basis of our algorithm 
to compute m(A). 

We choose a point X > such that for x E [X, oo) we may develop an asymptotic expansion 
for vfj(x,\), together with a precise estimate on the error committed. This expansion enables us 
to determine intervals in which ip{X, A) and if}'(X, A) lie, thus providing initial data to an interval 
based initial value solver that is used to compute complex intervals which enclose tp(a, A) and 
i/j'(a,\). The result ( |2. 3| ) yields an interval which encloses m(A). In section 3 we review the 
asymptotic method and interval ODE solver that is used to perform these tasks. 



3 Overview of the components of the algorithm 

3.1 Asymptotic theory 

The method that we use to obtain the asymptotic solution ( |1. 2| ) as x — > oo is the repeated 
diagonalization method of Eastham which is fully explained in the book || and here we shall 
be brief. The method is concerned with estimating and improving error terms in the asymptotic 
solution of the linear differential system 

Z'(x) = p(x){D + R{x)}Z{x) {a<x<oo) (3.1) 

where Z is an n— component vector, p is a real or complex scalar factor, D is a constant diagonal 
matrix 

D = dg(di,d 2 , -,dn) 
with distinct dk and R is a perturbation such that 

R(x) = 0(x~ s ) (x -> oo) (3. 2) 



for some 5 > 0. 

If it is the case that pR G L(a, oo), the Levinson asymptotic theorem can be applied to ( |3. 1| ) 
to give solutions 

rx 

Z k (x) = {e k + 7] k (x)} exp(d k / p(i)<&) (3.3) 

where e k is the unit coordinate vector in the k— direction and rj k (x) — > as x — > oo. The size of 
the error term is related to the size of i? as x — ► oo, and therefore the accuracy of ( 3. 3| ) can be 



improved if the perturbation R can be reduced as x — > oo. Under suitable conditions on p and 
R this improvement can be achieved by a sequence of repeated transformations which lead to a 



computational procedure to estimate the solution of (1. 2) together with a bound on the associated 
error. 

The sequence of transformations may be obtained either by an exact diagonalization or by an 
approximate diagonalization procedure. These methods are discussed in detail in ]lfl]. The exact 
diagonalization method involves the explicit construction of an n x n matrix T such that 

T-\x){D + R(x)}T(x) = D^x) 

and this in turn requires the explicit eigenvectors of D + R(x) which, although available for the 
second order system, are not generally known for the n— th order system. In this investigation we 
choose to work with the more generally applicable approximate diagonalization method which may 
be used for n— th order systems of differential equations. A discussion on the asymptotic method 



of exact diagonalization as applied to estimating the m function can be found in |10]. 



3.1.1 Approximate diagonalization 

We assume that R is a differentiate n x n matrix satisfying ( [3. 2j ), and we define an n x n matrix 
Pby 

PD - DP = R- dgR 
with diagonal entries pa = and other entries 

Pij = rij/(dj - di). (3. 4) 

We note that P = O(R) = 0(x~ s ) and the construction of the P matrix cancels out the dominant 
terms in the following system (|3. 5|) . With Z = (I + P)W, ( |3. 1| ) is transformed into the system 



W' = pip + S)W 



where we have written the n x n matrices 



D = D + dgR 

S={I + P)~ 1 (RP-PdgR-p~ 1 P') = 0(x- 2S ) +0(p- 1 x- s ~ 1 ) (3.6) 

which is of smaller magnitude than R. We write 

(/ + P)- 1 = I - p + P 2 + ... + {-\yp v + (-l) u+1 {I + P)- 1 ^^ 1 (3. 7) 

and specify an order of magnitude 0(x~ K ) which we wish to achieve as an error in the asymptotic 
solution of ( |3. 1| ). Substituting ( 3. 7 ) into ( 3, 6 ), we have 



S = V 2 + ...Vm-x + E (3.8) 

where 

V m = 0{x~ mS ) 
E = 0{x~ K ) 

and (M — 1)5 < K < MS. Thus v is chosen so that P u+1 in ( |3. 7[) gives rise to terms which 
contribute to E by at most this order of magnitude, and will be estimated in the final stages of the 
algorithm. 

This transformation procedure may be repeated for the W system ( [3. 5| ), but with a new P 
defined in terms of Vi which replaces R in Q3. We continue to use the matrix D and not 
D to simplify the construction of an efficient computational algorithm. However this procedure 
introduces additional terms into the analysis which must be eliminated at subsequent iterations of 
the algorithm. A repetition of the above process leads to a new matrix S viz. 

S = V 3 + ... + V M -i + E 

with new T^'s and a new E. 

The above ideas may be used to form the basis of an iterative procedure, which can be imple- 
mented in the symbolic algebra system Mathematical to compute the asymptotic solutions of (|3] 
~ip . Taking ( |3. 1[ ) as a starting point with m = 1, we have at the m— th stage 

— p(E m -\- R m )Z mi 
Rm — 'lm ^2m ■■■ *M—m,m E m , 

V jm = 0(x-( m+j -V s ), 



dgV lm = 0, 

D m = D + A m , 

Z m = (I + P m )Z m +i 

where P m is defined explicitly in terms of V\ m and D as in ( p. 4| ). Thus V\ m is eliminated at this 
stage. 

At the end of the process all the V's are eliminated and this gives 

Z'm = p(Dm + Rm)Zm 

where 

Rm = E M = 0(x- K ) 

and Mb > K. The Levinson theorem then yields the solution of the Zm system, and reversal of 
the M — 1 transformations gives 

z = {nf -\i + p m )}( e + v ) ex P ( / x p(...)dt) 

J a 

with 77 = 0(x~ K ), ( see [1C] for further details). 



3.2 Interval ODE solver 

In this sub-section we introduce briefly the concepts of interval arithmetic that we need to give a 
short account of Lohner's AWA algorithm. For an in-depth discussion of interval arithmetic, see 



[11 1, while Lohner's AWA algorithm is discussed in [12] and 



Denoting any of the four basic arithmetic operations by *, we define, for real intervals [a], [b], 

[a] ★ [b] = {a * b \ a e [a], b E [&]}. 

Thus we can compute an enclosure for [a]* [b] by obtaining computable upper, and lower bound, for 
[a] -k [b] which is derived from the lower and upper bounds of [a], [b] respectively, by some directed 
rounding facilities. Any algorithm that is realised on a computer consists of finitely many operations 
★ and thus an enclosure for the results of arithmetic operations which constitute the algorithm may 
be computed. In practice this simple approach would soon lead to an explosion of the interval 
width but many sophisticated techniques are available to control this phenomenon |TT] ]. 

Lohner's approach to computing an enclosure of the solution of initial value problems is based 
on the well known Taylor method for solving initial value problems. Suppose that a solution of the 



IVP 

u = f(x,u), u(0)=u , (3.9) 

where / : [0, oo) x BP — ► R n is sufficiently smooth, is known at some point xq. Then the solution 
at xq + h is 

u(x + h) = u{xq) + h(j)(x , h) + z X0+h (3. 10) 

where u(xo) +h(j)(xQ, h) is the (r — l)-th degree Taylor polynomial of u expanded about xq and z XQ+ h 
is the associated local error. This method lends itself well to computation since the coefficients 
of the polynomial may be computed via an automatic differentiation package by differentiating 



the differential equation ( 3. 9 ). However the error term is not known exactly since the standard 



formulae give, for some unknown r, 

z Xo+h = u<r\r)h r /rl, re[x ,x + h]. (3.11) 

In Lohner's algorithm, ( |3. 10j ) is used to advance an enclosure [m(xo)] for the solution u at xq, to 
one for the solution u at xq + h which we denote by [u(xq + h)]. A suitable enclosure for the error 
(Oil) is 

[z xo+h ]=f^([x ,x + h],[u])h r /rl 

provided that an enclosure [it] for {u(x) : xq < x < xq + h} can be computed. This is achieved by 
the following means. Choose some interval [u°] D [u(xo)] and try to prove that 

[u] = [u(x Q )] + [0,h]- f([x ,x + h],[u ]) C [A 

If this is true then Banach's fixed-point theorem implies that [u] is an enclosure for u(x) for all 
x £ (xo,xq + h). In order to achieve efficient performance and tight bounds, the details of the 
algorithm are more complex than this short overview can show. We refer the reader to JO] and |l]] 
for a complete discussion of the method. 



4 Results for q = —x a : a = 1, 2 

In this section we discuss the computation of m when a = and p = w = 1 and the potential 
q = —x a , < a < 2. In terms of the asymptotic analysis presented in section 3 this means 
that we take 5 = a/2. However, while the general algorithm is applicable to all a in this range, 
the implementation of Lohner's AWA interval ODE solver requires at least two derivatives of the 
function q, see ( |1. 1[ ), to be available at x = 0. Clearly this is not possible for < a < 2, a / 1, 



and for a in this range a revised algorithm is presented in section 5. Here we present an algorithm 
to compute m when q = —x or q = —x 2 . We remark that the algorithm which we present here is 
also applicable to problems where q = x a , 2 < a or a = 1, while that in section 5 covers the case 
< a < 2. 

We first write 2|) as the system 

1 
q-X 

and as in ]|, chapter 2] introduce the transformation 

/ 






This enables us to write Ql. 2j ) in the form ( |3. 1| ) with 



p = V^A, D = dg(l,-1) 

and 



R 



f-1 1 



4(q - A) 3 / 2 



V 



1 -1 



We next apply 6 iterations of the asymptotic algorithm of section 3.1 to obtain bounds on ip(x, A) 
and tp'(x,X) (X < x < oo) the L 2 [a, oo) solution obtained from the asymptotic algorithm. This 
gives intervals which enclose ip(X, A) and ip'(X, A) which are the initial data required by the AWA 
algorithm. The asymptotic analysis and estimation of the error is performed using purpose written 



Mathematica code, the detail of which is fully reported on in [10]. 

The C-XSC implementation of Lohner's algorithm is used with purpose written shell script to 
interface the asymptotic results to the interval arithmetic code. 

In all cases the enclosures that we obtain for m(A) are in agreement with the closed form of 
m(-) given in terms of either gamma functions or Bessel functions [^3|, and evaluated by numerical 
routines, (see [ |14| | and [|15| for further details). We further remark that the algorithm reported on 
in [pi] could not perform the computation required to produce the above results. 



A 


m(A) 


-1 + i 


n 7OQ7QQ7644403761 , n /|987r>7 119 0558884 : 
U.IZ.OI o-J 682 2 5801 2o T 1 U 1 0519233491 1 


i 


n wn5n 3090709121 -i- n fifi^-?fin 7324705136 i 

u.ooouou 2364882611 -h u.ooodou 6535858522 1 


0.5 + i 


n zL9zLru*?9 898149527 j- n 7nnozn 1148568185 ; 

U - 4z4U4 <J z 260330653 + U - < UUy41 0374761820 1 


0.1 + 0.1 i 


n ^000730075913471 , 1 nQfifiQ 4075678773 7 ; 

u - oyyy '29027001456 ' 1 - uyDDy 39419214239 1 


1 + i 


0-3 1 2 7 2 62 385936455 + 0-688666l34|g5 7 5 7 l i 


10 + 10 i 


n ini7^^ 71 39152363 1 n 9/l/l/l/lie474486887 : 
V.WU&6 6g73947S3g + U.z44441» 236960981 1 


1 + 0.5 i 


n 9RKQ793056163538 , n 7Q/L7C,'?8 875505129 ; 
U.ZDOy ( ^2659812543 + UJa4 ' OJ,5 307529690 1 


1 + 0.1 i 


n 1 S7Qc,«8032104862 , n 88 K7KA 9186094842 • 
u - lo ' OOD 7636023492 "+" u - 0Cl0 ' oq: 8512935614 1 


1 + 0.01 i 


0.1628908fl?!f!??I + 0.9047931°?™ i 


1 + 0.001 i 


0.1602956|l?gSI 7 i§ + 0.906631 ii 2 lii i 


1 + 0.0001 i 


0.160034658iia + 0.906814 ili 4 I« 4 i 


1 + 0.00001 i 


0.1600085f 447 o 24 i + 0.906832 4 l 22 i 7 |?I i 



Table 1: X = 10 and a = 1 , where e 6 (10) = 1.09576673x10 



A 


m(A) 


-1 + i 


0.72 1 54 63 224474 ^ 4 + 0.3676480? 4 ™ i 


i 


0.62 6 6 5 70 722 iHi 4 + 0.6266570 7 ™ 1 i 


0.5 + i 


0.4999306i |? 4 f 78 + 0.7052934^^94993 { 


0.1 + 0.1 i 


0.8975088111™ + 1.0328959K?i§S i 


1 + i 


0.367648°?—?? + 0.7215463?li 7 ||IfB i 


10 + 10 i 


0.10206645l 4 i 7 7 l? + 0.24554208™; 7 i 


1 + 0.5 i 


0.3372315g65?|4ggl + 0.86795375o55| 8 oo3 1 


1 + 0.1 i 


n 9 £«1 960389488290 , i 099891 6685241975 : 
u - ZODiz 59314609572 J - UZZOZ1 5104881830 1 


1 + 0.01 i 


n 99c;9798^ 5 0846718 , i nR1 nO"} 7010720598 i 
u - zzoz ' ZO °03204735 + 1 - UDluy,3 6911651063 1 


1 + 0.001 i 


0.22185279li 7 37 i + 1.0649430|™ 2 i 


1 + 0.0001 i 


n 221 ^n79 705977203 a. i nfiwsn 561152201 i 

U.ZZ1DU '^658584275 + 1 ■ UOOJZOU 461491257 1 


1 + 0.00001 i 


0.2214726f 7 42 i 28 7 § + 1.0653665il 474 flH i 



Table 2: X = 10 and a = 2 , where e 6 (10) = 1.20325893 x 10 



A 


m(A) 


-1 + i 


0.723783723^1§§| + 0.428707085!$^ i 


i 


n c;c;t;n505827338581 , n fifie,Qf:09723872275 : 
U - OOOU 499628151809 u - OOOJO 04136894624 1 


0.5 + i 


0.42404325^64205 + .70094107 5 7 S? i 


0.1 + 0.1 i 


n 6000855062924143 , i nQfi7992302211012 ; 
u -5998604127881420 ' i - UJD 5887988999492 1 


1 + i 


0.31272625?— + 0.688666159™ i 


10 + 10 i 


0.1017537056^4222 + 0.24444183557?^ i 


1 + 0.5 i 


n ORKQ734567852943 , n 7Q/1 7 c 53834479725 : 
U.ZDOy^ 1H49849902 + u - ' y4 ' °23351057398 1 


1 + 0.1 i 


1 «7qt;R8022703960 1 n oor;7E / |9051312494 : 
U.IO i oOD 765 3 04g066 -t- U.000 1 0^8667747971 1 


1 + 0.01 i 


n 1 «O«Qn9 12 0955266 , n 00/170^5903178042 : 

u.ib28yo 760961870 8 + u.yo4( , y33 4 84 6 3i 10 5 1 


1 + 0.001 i 


n ifin90 c ,fi 663626111 4- n onfifi^i 1260873454 ; 

U.10U/9t)D 5 g 74 3 4 i 46 + U.yUOOdli 165602 242 1 


1 + 0.0001 i 


O 1fiOO3/1 862442 0599 , n QOR«1 4 5096602290 • 
U - 1DUU,:54 4397312929 + u - yuD ° i 37578832769 1 


1 + 0.00001 i 


O 1 ROOO85817370011 , n OOK8395038096614 : 

U.lb00 08 4 88 253 4 979 + U.y0b832 3509gl0573 l 



Table 3: X = 40 and a = 1 , where e 6 (40) = 5.21570339 x 10" 15 . 

5 Results for q = — x a ,0 < a < 2, a = 0,p = w = 1 

As we mentioned at the beginning of section 4 we have to modify our algorithm to deal with values 
of a other than 1 and 2. We do this by choosing some number e > 0. Then, instead of solving 
2| ) over the whole of [0, X], we now solve the equation over the interval [e, X]. The following 
Theorem 5.1 enables us to obtain an enclosure for y(0) and y' (0) in terms of the enclosures for y(e) 
and y'(e). 

Theorem 5.1 Let c(x) = q(x) — A and, for some e > 0, let f,g 6 Ci[0, e] satisfy /(e) = y(e), 
g(e) = y (e). In addition let 

b = e [ | c | dt < 1. 
Jo 

Then 
1. 

I y(0) - /(0) |< jZTblf* \f'-9\dt + e£\g'-cf\dt] (5. 1) 

2. 

I y'(0) - g(0) \< —?-=-[ f \c\dt f | /' -g\dt+ f \ g - cf \ dt}. (5. 2) 
l — o Jo Jo Jo 



Proof: Define 



Then u(e) = and 



u{x) 



y{x) 

y'(x) 



u (x) 



y 0*0 

^ c(x)y(x) 

1 

c(x) 



/(*) 

g{x) 

u(x) — r(x), 



where 



Define 



f-9 
9 ~ cf 

T:(C[0, e ]) 2 -(C[0, e ]) 2 



by 



(Tv)(x) :-- 



1 



c(t) 

for x £ [0, e] and v £ (C[0, e]) 2 . Then it follows that 



)(t)dt- J r(t)dt, 



u = Tu. 



We shall prove that T has a globally unique fixed point in 



U :={v£ (C[0,e]) 2 :| vt(x) |< ai, | u 2 (x) \< a 2 (0 < x < e)} 



where aj and a 2 denote the right-hand sides of the inequalities ( |5. 1| ) and ( 5. 2| ), respectively. To 
show this we use the Banach fixed point theorem. This requires us to show 

1. T is a contraction with respect to some suitable norm; 

2. TU C U. 



In order to show (1) let 



v := max 



cg[0ie] {max{^ | Vi(x) |, 
V e 



Jo I c I dt 



V2(X) |}}. 



Since for all v,v € (C[0, e]) and x G [0, e], 



(Tv - Tv)i(x) | < y \v 2 {t)-v 2 {t)\dt<e\\v-v\\ ^ J \c\dt, 
(Tv-Tv) 2 (x) | < / | c(t) 1 1 vi (*) | dt< Veil 

u — £> || / | c | dt. 
•/.-• ' Jo 



we have 

|| Tv-Tv \\< ^le I \ c \ dt \\ v - v \\ . 

Since 



o 



Jo 



c I dt < 1 



it follows that T is a contraction. 

In order to establish part 2 we take v e U. Then 



(Tu)i(a;) | < ( \ v 2 \dt+ [ \ n | dt 

J X J X 

< ea 2 + I | r\ | dt 
= ai 



o 



and 

I (Tv) 2 (x) | < f | c || «i | dt+ [ \r 2 \dt 

J X J X 

< a\ | c | dt + / | T2 | (it 

JO JO 

= «2- 

Thus T has a globally unique fixed point in U which implies u £ U. In particular, 

I ^i(O) |< ai, | u 2 (0) |< a 2 

which establishes the theorem. 
5.1 Numerical examples 

We now turn to two examples which illustrate the use of Theorem 5.1 to compute enclosures for 
m. These examples are: p = w = 1, and 



1. q(x) = 

2. q(x) = -x 3 / 2 . 



In both these examples there is insufficient smoothness in the function q for Lohner's code to 
compute enclosures at x = 0. We therefore compute enclosures at x = e and use the above theorem 
to obtain enclosures at x = 0. 

In the first example we choose / and g to be the first-order Taylor polynomial approximations 
to y and y expanded about x = e. Writing Aq = y(e) and A\ = y '(e) these are respectively 

f(x) = A + (x-e)A 1 , 

g(x) = A 1 + (x-e)(-Ve-X)A . 



Since for this example 



[ e \ c \dt < -e 3 / 2 +|A|e 
Jo 3 



c - c(e) | dt = -e 3 / 2 
3 

c | (e-t)dt < -^e 5/2 + i | A | e 2 , 
15 2 

the inequalities (|5. 1| ) and ( 5. 2| ) become 



o 



y(0) - (A - eAt) \ < 



e 2 



1-6 2 (| A|+|^) 

{(^|A| + ^)Mo| + ^|A| + l^)|^|} 



2/ '(0)-(A 1 + e(A + v / ^)^o) I < 



F 3/2 



1-6 2 (| 

{[^ 3/2 (l A | +^)(| A | + 2 -V~e) + \] I A I +Ve(l \ A | +^) | A 1 |} 



In the example that we report on, A = 1 + i and we have taken e = 0.000015625. However, since 
Lohner's code uses a fixed step size algorithm, in order to achieve such a small value of e we have 
integrated over [10, 0.03125] with a step size of 0.03125, then reduced the step size to 0.000015625 
to perform the integration over [0.03125,0.000015625]. This yields an enclosure 

m (l +%) = 0.25936| 2 + 0.6719ffi 

We have compared our result with the estimate of m(l + i) obtained from the Brown Kir by Pryce 
code Q. That algorithm gives m(l + i) = 0.25937860 + 0.67196464i which is slightly outside our 
enclosures and reflects the lack of smoothness in q experienced by the Runga-Kutta method they 
employ. 



In the next example we take q = — x 3//2 and as before we take / and g to be the first-order 
Taylor approximations of y and y respectively, expanded about x = e. This time we get 

| y(0) - (A - eAi) | < 



y'(0) - (A 1 + e(A + e 3/2 )A ) \ < 



l-e 2 (| A | +§e 3 / 2 ) 

«i|A| + H e 3/ 2)Mo | +£( i m+ l £ 3/ 2)Ml | } 

r 2 



1 - e 2 (| A | +f e 3 / 2 ) 



{[^(1 A I + e 3/2 )(| A | +\e^) + ^ I A | +(I | A | +| e 3 / 2 ) | ^ |}. 



This gives an enclosure 



m(l + i) = 0.345229?!] + 0.705227^1 



We have investigated the possibility of choosing / and g to be the second-order Taylor polyno- 
mials. However, for both the examples that we have considered there appears to be no improvement 
in the bound over that achieved by linear functions /, g. 
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